﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using DotSpatial.Topology;
using DotSpatial.Data;

namespace GeostatisticalTool.Classes
{
	public class Utils
	{
       //double Pi=Math.PI; 
       //double  PiD2=Math.PI/2;
       //double PiD4=Math.PI/4; 
       //double Pi2=2*Math.PI;
       //double e=2.718281828459045235; 
       //double e10 = 1.105170918075647625;
       //double Deg=180/Math.PI;


       // public static double ChiSq(double x, double n) { 

       //     if(n==1 & x>1000) {return 0;}
       //     if(x>1000 | n>1000) { 

       //             var q=ChiSq((x-n)*(x-n)/(2*n),1)/2;

       //             if (x > n) { return q; } { return 1 - q; }
       //             }

       //     var Pi =Math.PI;
       //     var p=Math.Exp(-0.5*x); if((n%2)==1) { p=p*Math.Sqrt(2*x/Pi); }
       //     var k=n; while(k>=2) { p=p*x/k; k=k-2; }
       //     var t=p; var a=n; while(t>0.0000000001*p) { a=a+2; t=t*x/a; p=p+t; }

       //     return 1-p;

       // }
       // public static double Fmt(double y) { 
       // double  v ;
       // if(y>=0) { v=(y+0.0005); } else { v=(y-0.0005); }

       // return v;
       // }

        public static double ToRad(double value)
        {
            return value * Math.PI / 180;
        }
        public static double ToDeg(double value)
        {
            return value * 180 / Math.PI;
        }

       public static IFeatureSet FeatureFromPoints(List<Kpoint> points, IFeatureSet refe)
        {
            foreach(Kpoint p in points)
               {
                }

            return refe;
        }


        public static IGeometry DrawEllipse(Kpoint center, double majorRadius, double minorRadius, double azimuth)
        {
            int numVertex = 180;
            Coordinate[] coorPol = new Coordinate[numVertex+1];

            if (majorRadius == minorRadius)
            {
                int j = 0;
                for (int i = 0; i < 360; i += (360/numVertex))
                {
                    Kpoint v = Utils.AzimuthDist(center, i * Math.PI / 180.0, majorRadius);
                    coorPol[j] = new Coordinate(v.X, v.Y);
                    j++;
                }


                coorPol[numVertex] = coorPol[0];
            }
            else
            {
                int j = 0;
                for (int i = 0; i < 360; i += (360 / numVertex))
                {
                    double r = majorRadius * minorRadius / Math.Sqrt(
                        Math.Pow(minorRadius * Math.Cos(-1*azimuth + (i * Math.PI / 180.0)), 2) +
                         Math.Pow(majorRadius * Math.Sin(-1*azimuth + (i * Math.PI / 180.0)), 2)
                        );

                    Kpoint v = Utils.AzimuthDist(center, i * Math.PI / 180.0, r);
                    coorPol[j] = new Coordinate(v.X, v.Y);
                    j++;

                }
                coorPol[numVertex] = coorPol[0];

            }

            LinearRing pol = new LinearRing(coorPol);
            Polygon poly = new Polygon(pol);
            return poly;
        }


        public static List<IGeometry> DrawAxisEllipse(Kpoint center, double majorRadius, double minorRadius, double azimuth)
        {
            List<IGeometry> list = new List<IGeometry>();

            
            Coordinate[] east = new Coordinate[2];
            Coordinate[] north = new Coordinate[2];
            Coordinate centerP = new Coordinate(center.X,center.Y);


                east[0] = AzimuthDist(centerP, azimuth + Math.PI / 2, minorRadius);
                east[1] = AzimuthDist(centerP, azimuth - Math.PI / 2, minorRadius);
                north[0] = AzimuthDist(centerP, azimuth, majorRadius);
                north[1] = AzimuthDist(centerP, azimuth + Math.PI , majorRadius);


                list.Add(new LineString(east));
                list.Add(new LineString(north));
                return list;

        }




        public static double GetDistance(double x, double y, double z, Kpoint point)
        {

            return Math.Sqrt(((point.X - x) * (point.X - x)) +
                                 ((point.Y - y) * (point.Y - y)) + ((point.Z - z) * (point.Z - z)));
        }

        public static double GetDistanceXY(double x, double y,  Kpoint point)
        {

            return Math.Sqrt(((point.X - x) * (point.X - x)) +
                                 ((point.Y - y) * (point.Y - y)) );
        }

        public static double GetDistanceXY(Kpoint pointi, Kpoint point)
        {

            return Math.Sqrt(((point.X - pointi.X) * (point.X - pointi.X)) +
                                 ((point.Y - pointi.Y) * (point.Y - pointi.Y)));
        }

        public static double GetDistance(Kpoint pointi, Kpoint point)
        {

            return Math.Sqrt(((point.X - pointi.X) * (point.X - pointi.X)) +
                                 ((point.Y - pointi.Y) * (point.Y - pointi.Y)) + ((point.Z - pointi.Z) * (point.Z - pointi.Z)));
        }
        public static double GetSemivariance(Kpoint pointi, Kpoint point)
        {
            return 0.5 * (pointi.W - point.W) * (pointi.W - point.W);
        }

        public static double GetAzimuth(Kpoint center, Kpoint target)
        {
            return GetAzimuth(center.X, center.Y, target);
        }

        public static double[] GetAzimuth3D(Kpoint center, Kpoint target)
        {
            return GetAzimuth3D(center.X, center.Y, center.Z, target);
        }

        public static double GetAzimuth(double x, double y, Kpoint pi)
        {
            //double dx = (this.cpoints[j, 0] - this.cpoints[i, 0]);
            //double dy = (this.cpoints[j, 1] - this.cpoints[i, 1]);
            double dx = pi.X - x;
            double dy = pi.Y - y;
            double az = double.NaN;
            if (dx == 0 && dy > 0) az= 0;
            if (dx == 0 && dy < 0) az= (Math.PI);
            if (dx > 0 && dy == 0) az= (Math.PI / 2);
            if (dx < 0 && dy == 0) az= (3 * Math.PI / 2);
            if (dx > 0 && dy > 0) az= (Math.Atan(dx / dy));
            if (dx > 0 && dy < 0) az= (Math.PI + Math.Atan(dx / dy));
            if (dx < 0 && dy < 0) az= ((Math.PI) + Math.Atan(dx / dy));
            if (dx < 0 && dy > 0) az= ((2 * Math.PI) + Math.Atan(dx / dy));

            if (az > (Math.PI * 2))
            {
                double val = (az) / (Math.PI * 2);
                return (val - Math.Floor(val)) * (Math.PI * 2);
            }

            return az;
        }

        public static double[] GetAzimuth3D(double x, double y, double z, Kpoint pi)
        {

            double azXY = GetAzimuth(x, y, pi);
            double dx = GetDistanceXY(x, y, pi);
            double dy = pi.Z - z;
            double az = double.NaN;
            if (dx == 0 && dy > 0) az = 0;
            if (dx == 0 && dy < 0) az = (Math.PI);
            if (dx > 0 && dy == 0) az = (Math.PI / 2);
            if (dx < 0 && dy == 0) az = (3 * Math.PI / 2);
            if (dx > 0 && dy > 0) az = (Math.Atan(dx / dy));
            if (dx > 0 && dy < 0) az = (Math.PI + Math.Atan(dx / dy));
            if (dx < 0 && dy < 0) az = ((Math.PI) + Math.Atan(dx / dy));
            if (dx < 0 && dy > 0) az = ((2 * Math.PI) + Math.Atan(dx / dy));

            if (az > (Math.PI * 2))
            {
                double val = (az) / (Math.PI * 2);
                az= (val - Math.Floor(val)) * (Math.PI * 2);
            }

            return new double[2]{azXY, az};
            
        }



        public static double GetAzimut(Coordinate origin, Coordinate target)
        {
            double dx = target.X - origin.X;
            double dy = target.Y - origin.Y;

            // if (dx == 0 && dy == 0) return 0;
            if (dx == 0 && dy > 0) return 0;
            if (dx == 0 && dy < 0) return Math.PI;
            if (dx > 0 && dy == 0) return Math.PI / 2;
            if (dx < 0 && dy == 0) return 3 * Math.PI / 2;
            if (dx > 0 && dy > 0) return Math.Atan(dx / dy);
            if (dx > 0 && dy < 0) return Math.PI + Math.Atan(dx / dy);
            if (dx < 0 && dy < 0) return (Math.PI) + Math.Atan(dx / dy);
            if (dx < 0 && dy > 0) return (2 * Math.PI) + Math.Atan(dx / dy);
            return 0;
            //  return Math.Atan(dx / dy);

        }

        public static double GetAzimut(double x, double y, double x1, double y1)
        {
            double dx = x1 - x;
            double dy = y1 - y;

            // if (dx == 0 && dy == 0) return 0;
            if (dx == 0 && dy > 0) return 0;
            if (dx == 0 && dy < 0) return Math.PI;
            if (dx > 0 && dy == 0) return Math.PI / 2;
            if (dx < 0 && dy == 0) return 3 * Math.PI / 2;
            if (dx > 0 && dy > 0) return Math.Atan(dx / dy);
            if (dx > 0 && dy < 0) return Math.PI + Math.Atan(dx / dy);
            if (dx < 0 && dy < 0) return (Math.PI) + Math.Atan(dx / dy);
            if (dx < 0 && dy > 0) return (2 * Math.PI) + Math.Atan(dx / dy);
            return 0;
            //  return Math.Atan(dx / dy);

        }

        public static double VerifyAzimuth(double az)
       {
           if (az > (Math.PI * 2))
           {
               double val = (az) / (Math.PI * 2);
               return (val - Math.Floor(val)) * (Math.PI * 2);
           }
           return az;
        
        }
        
        
        public static double SumAzimuth(double value, double value2)
        {

            return VerifyAzimuth(value + value2);

        }

        public static Kpoint AzimuthDist(double x, double y, double azi, double dist)
        {
            return new Kpoint(x + Math.Sin(azi) * dist, y + Math.Cos(azi) * dist, 0, 0);
        }

        public static Coordinate AzimuthDistC(double x, double y, double azi, double dist)
        {
            return new Coordinate(x + Math.Sin(azi) * dist, y + Math.Cos(azi) * dist);
        }

        public static Kpoint AzimuthDist(Kpoint ptoI, double azi, double dist)
        {
            return new Kpoint(ptoI.X + Math.Sin(azi) * dist, ptoI.Y + Math.Cos(azi) * dist, 0, 0);
        }

        public static Kpoint AzimuthDist(Kpoint ptoI, double[] azi, double dist)
        {
            double r = Math.Sin(azi[1]) * dist;

            return new Kpoint(ptoI.X + Math.Sin(azi[0]) * r, ptoI.Y + Math.Cos(azi[0]) * r, Math.Cos(azi[1]) * dist, 0);
        }



        public static Coordinate AzimuthDist(Coordinate ptoI, double azi, double dist)
        {
            return new Coordinate(ptoI.X + Math.Sin(azi) * dist, ptoI.Y + Math.Cos(azi) * dist);
        }

		}

	public struct Stat
	{

		public double value;
		public long num;
		public double quad;
		private double max;
		private double min;
		private long posMin;
		private long posMax;

		public long PosMin
		{
			get
			{
				return this.posMin;
			}
		}

		public long PosMax
		{
			get
			{
				return this.posMax;
			}
		}


		public double Min
		{
			get
			{
				return this.min;
			}
		}


		public double Max
		{
			get
			{
				return this.max;
			}
		}

		public double Mean
		{

			get
			{
				if (num > 0)
					return value / this.num;
				else
					return 0;
			}
		}

		public double StdDevP
		{

			get
			{
				if (num > 0)
					return Math.Sqrt((quad - (2 * this.value * this.Mean) + (this.num * this.Mean * this.Mean)) / this.num);
				else
					return 0;
			}
		}

		public double StdDevS
		{

			get
			{
				if (num > 2)
					return Math.Sqrt((quad - (2 * this.value * this.Mean) + (this.num * this.Mean * this.Mean)) / (this.num - 1));
				else
					return 0;
			}
		}

		public Stat(bool v)
		{
			this.value = 0;
			this.quad = 0;
			this.min = double.MaxValue;
			this.max = double.MinValue;
			this.num = 0;
			this.posMax = 0;
			this.posMin = 0;
		}


		public Stat(double value)
		{
			this.value = value;
			this.quad = value * value;
			this.min = value;
			this.max = value;
			this.num = 1;
			this.posMax = 1;
			this.posMin = 1;


		}

		//Stat(double value, int v)
		//{
		//    this.value = value;
		//    this.quad = (value * value);
		//    this.num = v;
		//}

		Stat(Stat c, double value1, long v)
		{
			this.value = c.value + value1;
			this.quad = c.quad;
			this.quad += value1 * value1;
			this.num = v;
			//this.posMin =v;
			//this.posMax=v;


			if (c.Min > value1)
			{
				this.min = value1;
				this.posMin = v;
			}
			else
			{
				this.min = c.Min;
				this.posMin = c.posMin;
			}

			if (c.Max < value1)
			{
				this.max = value1;
				this.posMax = v;
			}
			else
			{
				this.max = c.Max;
				this.posMax = c.posMax;
			}

		}

		public static Stat operator +(Stat c1, Stat c2)
		{

			return new Stat(c1, c2.value, c1.num + 1);
		}

	}
	public struct Covariance
	{
		Stat vi;
		Stat vj;
		double value;

		public Covariance(bool v)
		{

			this.vi = new Stat(false);
			this.vj = new Stat(false);
			this.value = 0;
		}

		public Covariance(double vvi, double vvj)
		{
			this.value = vvi * vvj;
			this.vi = new Stat(false);
			this.vj = new Stat(false);
			this.vi += new Stat(vvi);
			this.vj += new Stat(vvj);
		}

		Covariance(Covariance cc1, double v1, Covariance cc2, double v2)
		{
			this.vi = cc1.vi + cc2.vi;
			this.vj = cc1.vj + cc2.vj;
			// this.value = cc1.value;
			this.value = v1 + v2;
		}

		public static Covariance operator +(Covariance c1, Covariance c2)
		{

			return new Covariance(c1, c1.value, c2, c2.value);
		}


		public double Result()
		{
			return this.value - this.vi.Mean * this.vj.Mean * this.vi.num;

		}
	}
    public struct Kpoint:ICloneable
    {
        public double X;
        public double Y;
        public double Z;
        public double W;

        public double Zinte;
        public double Error;
        public double SdtError;

        public long id;

        public Kpoint(double X, double Y, double Z,double W )
        {
            this.X = X;
            this.Y = Y;
            this.Z = Z;
            this.W = W;
            this.id = 0;
            this.Zinte = 0;
            this.Error = 0;
            this.SdtError = 0;
        }

        public Kpoint(double X, double Y, double Z, double W, long id)
        {
            this.X = X;
            this.Y = Y;
            this.W = W;
            this.Z = Z;
            this.id = id;
            this.Zinte = 0;
            this.Error = 0;
            this.SdtError = 0;
        }

        public object Clone()
        {
            Kpoint v= new Kpoint(this.X,this.Y,this.Z,this.W,this.id);
            v.Zinte = this.Zinte;
            v.Error = this.Error;
            v.SdtError = this.SdtError;
            return v;
        }

        public static bool operator ==(Kpoint c1, Kpoint c2)
        {

            return (c1.X == c2.X && c1.Y == c2.Y && c1.Z== c2.Z);
        }

        public static bool operator !=(Kpoint c1, Kpoint c2)
        {

            return !(c1== c2);
        }

        public override bool Equals(object obj)
        {
            if (obj.GetType() != typeof(Kpoint))
                return false;

            Kpoint p = (Kpoint)obj;
            return ((this.X == p.X) && (this.Y == p.Y) && (this.Z== p.Z));
        }

        public override int GetHashCode()
        {
            
            return Convert.ToInt32(X) ^ Convert.ToInt32(Y)^Convert.ToInt32(Z) ;
        }

    }

    public struct KDistance : IComparable
    {
        private double distance;

        public double Distance
        {
            get { return distance; }
            set { distance = value; }
        }
        public double azimuth;
        public Kpoint point;
        public KDistance(double distance, double azimuth, Kpoint p)
        {
            this.distance = distance;
            this.azimuth = azimuth;
            this.point = p;
        }

        public int CompareTo(object obj)
        {
            if (obj == null) return 1;
            KDistance otherDist = (KDistance)obj;
            if (otherDist != null)
                return this.distance.CompareTo(otherDist.distance);
            else
                throw new ArgumentException("Object is not a DistanceOrder");
        }

        public static bool operator ==(KDistance c1, KDistance c2)
        {

            return (c1.distance == c2.distance);
        }

        public static bool operator !=(KDistance c1, KDistance c2)
        {

            return !(c1 == c2);
        }
    }

    public struct RGB
    {
        public Stat R;
        public Stat B;
        public Stat G;

        public RGB(bool v)
        {
            this.R = new Stat(false);
            this.B = new Stat(false);
            this.G = new Stat(false);
        }

        public RGB(int r, int g, int b)
        {
            this.R = new Stat(false);
            this.B = new Stat(false);
            this.G = new Stat(false);
            this.R = this.R + new Stat(r);
            this.B = this.B + new Stat(b);
            this.G = this.G + new Stat(g);
        }

        RGB(RGB c, RGB d)
        {
            this.R = c.R + d.R;
            this.B = c.B + d.B;
            this.G = c.G + d.G;
        }

        public static RGB operator +(RGB c1, RGB c2)
        {

            return new RGB(c1, c2);
        }

    }

}
